###############################
# Analysis: Ruler survival
#
# Replication Material for:
# Continuity or Change? (In)direct Rule in British and French Colonial Africa
# 
# Carl Mueller-Crepon, 2020
# International Organization
#
# File Description:
# Analysis of ruler survival.
#
# Called from scripts/polities/analysis.R 
#
##################################

# Prepare Data
ruler.yrs.df$ruler.death <- ifelse(ruler.yrs.df$year == ruler.yrs.df$end.year & ruler.yrs.df$end.year!= 2006,
                                   1,0)
ruler.yrs.df$time <- ruler.yrs.df$year - ruler.yrs.df$start.year 
ruler.yrs.df$start <- ruler.yrs.df$time - 1
ruler.yrs.df$stop <- ruler.yrs.df$time
ruler.yrs.df$polity.age <- ruler.yrs.df$year - ruler.yrs.df$start.min
ruler.yrs.df$keep <- ifelse(ruler.yrs.df$end.year - ruler.yrs.df$start.year > 80 | grepl("unknown|colonial|ruled|direct|vacant", ruler.yrs.df$ruler) |
                              grepl("colonial", ruler.yrs.df$type) | ruler.yrs.df$year > 1962 |
                              !ruler.yrs.df$type.polity %in% c("",NA) |
                              (ruler.yrs.df$time == 0 & ruler.yrs.df$ruler.death == 1) | ruler.yrs.df$time < 1 | 
                              is.na(ruler.yrs.df$end.year) | is.na(ruler.yrs.df$start.year),
                            F, T)

# Get vars from polities
ruler.yrs.df <- join(ruler.yrs.df, polity.yrs.df[,c("polity.id","year",colnames(polity.yrs.df)[!colnames(polity.yrs.df) %in% colnames(ruler.yrs.df)])],
                     type = "left")

# Running variable to/from point of colonization
ruler.yrs.df$run.var <- ifelse(ruler.yrs.df$notcolonized == 0, ruler.yrs.df$start.col, ruler.yrs.df$yrs.2.col)

##########################
# Model setup

# ... base specification rulers
base.controls.rulers <- c( "col.brit","col.other", "notcolonized", "polity.age")

# ... variables to keep and labels
keep.controls.rulers <- c(1:length(base.controls.rulers))
keep.labels.rulers <- c("British rule","Other colonizer", "Not (yet) colonized","Polity age (log)")

# Strata
strata.ruler <- "strata(capital.id)"

##############################################
# Descriptive Plot ###########################
# Appendix Figure A10
##############################################

stub <- "descr.plot"

# Make plot data
plot.df <- ruler.yrs.df[ruler.yrs.df$year > 1500 & ruler.yrs.df$keep & ruler.yrs.df$action.type != "decolonization" &
                          ruler.yrs.df$polity.id %in% unique(polity.yrs.df$polity.id[polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1]),]
plot.df$groups <- ifelse(plot.df$notcolonized == 1, "Not colonized",
                         ifelse(plot.df$col.brit, "British rule",
                                ifelse(plot.df$col.frnc, "French rule", "Other colonizer")))
plot.df$run.var[plot.df$notcolonized == 1] <- plot.df$run.var[plot.df$notcolonized == 1] * -1
plot.df$Status <- plot.df$groups

# Plot survival against running variable
png(paste0(fig.path, "descr.rulers.run.png"), width = 7, height = 4, unit = "in", res = 300)
g <- ggplot(plot.df[plot.df$run.var > -50 & plot.df$run.var < 50,], 
       aes(x = run.var, y = ruler.death, group = Status, color = Status, pch = Status, lty = Status)) + 
  geom_smooth(se = F, method = "loess") +
  stat_summary_bin(fun.y='mean', bins=100, size=1.5, geom='point') +
  xlab("Time to colonial conquest") + ylab("Proportion of death/deposition of rulers")
print(g)
dev.off()

##################################
# Rulers: Main model #############
# Appendix Table A14 #############
##################################

stub <- "ruler.strata"

# ... running variables
run.vars <- c("run.var*notcolonized", "run.var*col.brit","run.var*col.other", "run.var*col.frnc")
run2.vars <- c("run.var*notcolonized", "run.var*col.brit","run.var*col.other", "run.var*col.frnc",
               "I(run.var^2)*notcolonized", "I(run.var^2)*col.brit","I(run.var^2)*col.other", "I(run.var^2)*col.frnc")

# ... full specification 
spec <- paste0("Surv(start, stop, ruler.death)", " ~  ", 
               c(paste(c( base.controls.rulers,strata.ruler), collapse = "+"),
                 paste(c(base.controls.rulers,
                         paste("I(", run.vars, ")"),
                         strata.ruler), collapse = "+"),
                 paste(c(base.controls.rulers,
                         paste("I(", run2.vars, ")"),
                         strata.ruler), collapse = "+"),
                 paste(c(base.controls.rulers, "polity.death",
                         paste("I(", run2.vars, ")"),
                         strata.ruler), collapse = "+")),
               "+ cluster(ruler.id)")



# Estimation
this.m <- lapply(spec, function(form.str){coxph( as.formula(form.str), 
                                                 ruler.yrs.df[ruler.yrs.df$year > 1500 & ruler.yrs.df$keep & ruler.yrs.df$action.type != "decolonization" &
                                                                ruler.yrs.df$polity.id %in% unique(polity.yrs.df$polity.id[polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1]),])})


# Stargazer

# ... Info
add.lines <- list(latex.any.row("Strata: ",c("capital", "capital", "capital", "capital") ),
                  latex.any.row("Running linear: ",c("no", "yes", "yes", "yes") ),
                  latex.any.row("Running quadratic: ",c("no", "no", "yes", "yes") ),
                  latex.any.row("Sample: ",c("post-1500", "post-1500", "post-1500", "post-1500") ))

fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(this.m, se = lapply(this.m, function(m){diag(m$var)^.5}),
            title="Death/deposition of rulers before and during colonial rule (1500--): Cox Proportional Hazards",
            dep.var.caption = "Death/deposition of ruler",dep.var.labels.include = FALSE,
            keep= c(keep.controls.rulers, 5), 
            covariate.labels = c(keep.labels.rulers,"End of line of succession"),
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,
            notes = nice_notes(c(notes.p,"Standard errors are clustered on the ruler-level.")), 
            notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)
